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Abstract 

We show that, when handled properly, the Cauchy principal value 
integral j^^^ f{x){x — T)~^dx {a < t < b) can be computed very easily 
and accurately using any reliable adaptive quadrature. We give de- 
tailed estimations of the roundoff errors in the case the function / has 
bounded first derivative in the interval [a, 6]. 
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1 Introduction 

We consider the problem of numerical evaluation of the Cauchy principal 
value integral 

IrMf) = f^dx = hm ( r\ f )I^dx, (1.1) 

Ja X-T 1,^0+ V J a Jr+fi J X-T 

where r G (a, b), and the function / has bounded first derivative. In general, 
the integral (II. ip exists if / is Holder continuous. The integrals of the 
type (jl.ip appear in many practical problems related to aerodynamics, wave 
propagation or fluid and fracture mechanics, mostly with relation to solving 
singular integral equations. 

A great many papers related to numerical evaluation of the integrals of 
the form p.ip have been published so far. Some of them are 121 [3l HI [6l [71 [8] . 
A nice survey on the subject, along with a large number of references, is 
presented in [H §2.12.8]. 

Even though so many algorithms have been proposed, the subroutines 
for computing the integrals of the type (II. 1|) are not commonly available in 
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systems for scientific computations. This is probably because most of the 
methods assume some properties of the function /, e.g. that / is analytic, 
/ is not rapidly oscillating etc. On the other hand, almost every system is 
equipped with one or more subroutines for automatic computation of the 
integrals 

[ f{x)dx. (1.2) 

J a 

These algorithms are usually based on the so called adaptive quadratures 
and they can compute the integrals of the form (|1.2p for a very wide range 
of integrands. The natural question arises: can an adaptive quadrature be 
used for computing the integral (jl.ip ? 

In this paper, we give the positive answer to that question. The next 
section contains the formulation of our algorithm. We apply two known an- 
alytical transformations that convert the integral (jl.ip into the sum of two 
nonsingular integrals. In Section [31 we discuss the problem of loss of signifi- 
cant digits that theoretically may appear when computing the transformed 
integral. In Section HI we present some numerical examples to validate the 
usefulness of the proposed method. 



2 Analitycal tranformations 

Without the loss of generality, we may restrict our attention to the case 
a = — 1 and h = 1. The computation of the Cauchy principal value integral 

IrU) ^ Ir,-lM) = f ^dx (2.1) 
J-lX-T 

may at first seem quite easy, if we observe that by a simple change of vari- 
ables, for 5 = min{l + r, 1 — r}, we obtain 



IrU) = / ^dx+ I 

J\x-t\>5 X — T 



x — t\>S X T 



5 X-T 



where we use the convention that f^^_^^yg = f^^s^ if = 1+t, and J[^_^[>5 = 

The formula (|2.2p was applied for the first time by Longman in [1], 
and it was derived by splitting the function / into the odd and even parts. 
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Both integrals on the right hand side of ()2.2p exist in the Riemann sense. 
We should note that if the function / has bounded first derivative in the 
neighbourhood of r, then this integral is not even singular. 

The first integral on the right hand side of (|2.2p was commonly ignored, 
as it is always a proper one. However, if r is close to —1 or 1, then this 
integral is near singular and standard quadratures may fail when applied 
directly. 

In many algorithms, another transformation of the integral Irif) is used, 
usually being called subtracting out the singularity. We have 



Irif) = / ^dx+ [ 



X — T 



(2.3) 



1 + T J_i X — T 



A direct application of the above formula is commonly not recommended 
(see, e.g., O |5]) due to possible severe cancellation, if a quadrature node 
happens to be very close to r. 

The author, however, has never seen the two above approaches being 
put together. Prom (12. 2p and (12. 3p we immediately obtain 

Irif) = fir) log[^+ / gix)dx+ [\ix)dx, (2.4) 

-L + T J|a;-r|><5 Jo 

where 

fix) - fir) /(r + x)-/(r-x) 

q(x) = and nix) = . 

^ ^ x-T ^ ' X 

If the function / has bounded first derivative, none of the integrals on the 
right hand side of (12. 4[) is singular or near singular (unless / itself has sin- 
gularities just outside the interval [—1,1]), and, when approximating these 
integrals numerically, the distance between r and any of the quadrature 
nodes is never smaller than b. 



3 Estimating roundoff errors 

In this section, we will show that the formula (j2.4p is numerically safe, i.e., 
we will prove that if the integral (j2.ip is approximated numerically using the 
equality (j2.4p . then the absolute error of the approximation, which results 
from the roundoff errors, is practically very small. 
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Let us denote by e the precision of the arithmetic used, and by (7e(x) 
and hi,{x) the numerically computed values of the functions g and h at the 
point X. For a moment, we will assume that the values of the function / are 
computed exactly. 

Lemma 1. If we define Di := max|2,|<]^ then for all \x\ < 6 

\he{x) - h{x)\ ^ + AeDi. (3.1) 

X 

Proof. We have 

, /((r(l + 7i) + x(l + 72))(l+73)) - /((t(1+7i) - x(l+72))(l + 74)) 
'^^^ x(l + 75) 

(3.2) 

where |7i|, I72I, I73I, I74I < e and I75I <, 2e. Further, 

(t(1+7i) + x(l + 72))(l + 73) ^ r + X + (71 + 73)t + (72 + 73)2; 

= r + X+ (r + x)/3 

for some |/3| ^ 2e, and, consequently, 

|/((r(l+7i) + x(l+72))(l+73))-/(r + x)| ^ \{t + x)/3\Di < 2eDi, (3.3) 

as |t + x| < 1. Analogously, 

|/((t(1+7i) - x(l+72))(l + 74)) - /(r - x)\ ^ 2eD^. (3.4) 

Now, from ([32]), and ([331) we obtain 

(/(r + x)-/(r-x))(l + 75) 



\he{x)-h{x)\ 



he{x) 
he{x) 



x(l + 75) 

{f{T + x)-f{T-x)) 75 (/(r + x)-/(r-x)) 



x(l + 75) 1 + 75 



□ 

Lemma 2. /f |a; — r| < 8e, then 

\g,{x)-g{x)\ ^ (3.5) 

\x — r 
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Proof. Similarly as in ()3.2p . we have 

/(x(1 + 7i))-/(t(1+72)) 



9e[X) 



(2;(l + 7i) -r(l+72))(l + 73)' 
where |7i|, I72I < e and I73I ^ 2e. Moreover, 

(x(l+7i) -r(l + 72))(l+73) - (x-r) (3 



(3.6) 



X — T 



(3 



X — T 



X — T 

for some \/3\ ^ 3e, which implies that 

(x(l + 7i) -t(1 + 72))(1 + 73) = {x-t)(i 
In addition, as |x|, |t| < 1, we have 

|/(x(l + 7i))-/(x)| <eDi and |/(r(l + 72))-/(r)| <eDi, 
which together with (j3.6j) and \(3\ ^ 3e, jx — r| < 8e gives 

fix) -fir) , P 



\ge{x) - g{x)\ 



Qeix) 



1 + 



X — T 



^ 2eDi + pPi ^ 8eDi 



\x — r 



□ 



In both lemmas we assumed that the values of the function / are com- 
puted exactly. In practice, this is obviously not true. However, if the com- 
putation of fix) for every x E [—1, 1] is numerically backward stable, then 
all the above results remain true, only the constant factors change. 

3.1 The first approach — cutting of the singularity 

From ()3.ip and (13. 5p we immediately obtain that 



Eeif) := [ h 

J\x~t\>S 



(x) — g{x)\dx + / \hs{x) — h{x)\dx 



rdx < IGeDi / —dx, 

-1 \x -t\ Jo x 



(3.7) 
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which may not look hke a good approximation of the cumulated roundoff 
errors. Thus, it seems quite natural to replace the last integral in ()2.4p by 

h{x) dx 



for some very small value of (/U < 5). In such a case, we have 

E,{f) ^ 16eI)ilog(/i-i). 

By A(f, a, b) we will denote the value of the integral (jl.2p approximated 
numerically by some algorithm A. We will also define 



EAif,a,b) :-- 



b 



f{x)dx - A{f,a,b) 



Now, as \ Jq h[x)dx\ < 2fiDi, we obtain 

)+A{ge,T + 6,l)+A{he,^l,6] 

(3.8) 



Irif) = f{r) \oE\-^ + A{g,,-l,T-5)+A{ge,T + 5,l) + A{he,^i,6)+£, 
i + r 



where 



1^1 <, EA{ge,-l,T-5)+EA{ge,T + 5,l)+EA{h„fI,5) 

+ 16eZ?ilog(/u-i) + 2^Z?i. 
3.2 The second approach — open-type quadratures 



If, for approximating the integrals on the right hand side of ()2.4p . we use 
a quadrature of the open type, then we are guaranteed that does not belong 
to the set of nodes of the quadrature A{h,0,5). Observe that, instead of 
(j3.8p . we may write 



Irif) = f{T) log^-^ + A{g,-1,T - 6) + A{g,T + 6,1) + A{h,0,6) + £2, 

i + T 

where 

1^2! ^ EA(,g,-l,T - 6) + EA{g,T + 6,1) + EA{h,0,6) 

+ A{\ge - g\,-l,r - 6) + A{\ge - g\,r + 6,1) + A{\he - h\,0,6). 

Analogously to ()3.7p . we have 

A{\g,-g\,-l,r-5)+A{\ge-g\,r+6,l)+A{\h,-h\,0,5) ^ 16eZ5iA(^, 0, 1). 
Now, we use the following 
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Observation 1. Assume that 

n-l 



1 i -j^ 

= Y,G{-,Sk,Sk+i) 

where 



X 



G{f, Sk, Sk+i) = ^ Bkjfixkj) 

j=0 

is a Gauss or Gauss-Kronrod quadrature rule applied to the interval [sk, Sfc+i]- 
If xofi denotes the smallest node of the quadrature rule G{f,so,si), then 

^(i,0,l) < C^„log(x^i), (3.9) 

and Cmn < 2 for all m > 1, n> 0. 

The above observation has not yet been proved. However, it was verified 
experimentally for 1 < m < 100 and thousands different values of n and 
sq,si, . . . , Sn- The greater the value of n is, the smaller is the constant Cmn- 
E.g., if m > 14, then Cmn < 1-3. 

From (j3.9p . we immediately obtain 

\£2\ ^ EAig, -1, T-6)+ Ea{9, t + 6,1)+ EAih, 0, 6) 
+ 32eDilog(xoi). 



4 Numerical experiments 
4.1 The algorithm 

Before we formulate our algorithm, we have to deal with one more problem, 
which, surprisingly, is usually ignored in the literature. Suppose that I7I < e. 
In some cases, the exact value of the integral may be significantly 

different from Irif)- It means that even if we would have a perfect algorithm 
for computing the integrals of the type (|2.ip , we may get wrong results only 
by rounding the parameter r. Such a situation takes place if |r| is close to 
1 or if /' changes rapidly in the neighbourhood of r. 

We can see that if |r| ~ 1, then small changes of r affect mostly the first 
term on the right hand side of ()2.4p . If we denote L{x) := log((l— x)/(l+3;)), 
then 

|L(r(l + 7))-^r)| ^ |rL'(r)| 



7 



and, consequently, 



t/(t)| 



(4.1) 



min{l + r, 1 — r} 



The error that results from the rapid changes of /' in the neighbourhood of 
r seems to be much more difficult to estimate theoretically. However, we 
have confirmed experimentally that it can be approximated as follows: 



for a properly selected constant C. 

The proposed algorithm was programmed in the Matlab language. As 
the algorithm A we used the build-in Matlab adaptive integrator quadgk. 
The approximation to the integral (j2.ip is computed as follows: 



In order to estimate the approximation error, we make use of the formula 
p.lOp . where, instead of EA{g, •, •) and EA{h, •, •), we use the error bounds 
for EAiQerr) and EA{hi;,-,-) provided, in fact, by the quadgk subroutine. 
Also, because the last term in (j3.10p is rather pessimistic and because we 
do not know the values of xo,o and we replaced the factor 32Di log(xQQ) 
by 8|/'(t)|. The final error estimation is increased by the formulas given in 
()4.ip and (|4.2p with C := 8. The values of /'(r) and /"(r) are approximated 
by simple divided differences. 

4.2 Numerical results 

As a counterpart we have chosen the algorithm presented in [3], which is 
based on the idea of Clenshaw-Curtis quadrature, and which we have verified 
to be fast and accurate. This algorithm is also of an automatic type, i.e., it 
tries to approximate the integral (12. ip within the prescribed error tolerance 
and also, together with the computed approximation, returns the estimated 
error bound. 

We performed our test for the following set of examples: 



(4.2) 



Irif ) ~ A{g, -1, r -6)+ A{g, t + 6,1)+ A{h, 0, 6). 



fix) 



T = 0.5 



fix) 



sin(5503;) 



r = 0.8 



fix) 



V^2 + cos(200x), T = 0.7, 
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f{x) = log2(1.0001 -x), T = 0.99, 

f{x) = V|cos(44x)|3, r = -0.6, 

f{x) = Vl -x2cos(100x), T = 0.5, 

f[x) = e^', T = 0.9999999, 

/(x) = 6-100(^+0-4)^ sin(e-io^), r = -0.41. 

By CC we denote the algorithm of [3j, while by A the algorithm based on 
the adaptive quadrature, proposed in this paper. In Matlab, e 1.1 • lO^i^, 
while the requested error tolerance was equal to 10 -i^. The experiments 
were performed on the computer with Intel Core 15 3.66 GHz processor. 
In tables, we give the values of the absolute errors of the computed ap- 
proximations, the estimated error bounds for these errors provided by the 
algorithms, and the computation times. 

Table 1: Comparison (absolute error, error estimation and computation time) of 
the algorithm of ^ (CC) and the present method (A) in the case of the integral 
J^-^e''{x-0.5)-^dx (IT = 0.0003s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


5.6 • 10-16 


1.4-10-15 


LOT 


A 


7.8 • 10-16 


3.8 • 10-15 


5.7r 



Table 2: Comparison (absolute error, error estimation and computation time) of 
the algorithm of [3] (CC) and the present method (A) in the case of the integral 
/^^ sin(550a;)(2: - 0.8)-^dx (IT = 0.001s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


1.9 • 10-14 


2.9 • 10-14 


LOT 


A 


2.8 • 10-13 


5.1 • 10-13 


4.7T 



In Tables [UH] we present the results in the case the function / is analytic 
in a complex region containing the interval [—1,1]. As the algorithm of 
[3] was designed for this class of functions, it performs a little faster then 
the algorithm presented in this paper. It is not a surprise, as adaptive 
quadratures are not meant to be daemons of speed, but to compute integrals 
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Table 3: Comparison (absolute error, error estimation and computation time) of 
the algorithm of [3] (CC) and the present method (A) in the case of the integral 
ill ^2 + cos(200a:)(.T - 0.7)-^dx (IT = 0.0029s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


5.3 • 10-15 


1.5-10-14 


LOT 


A 


2.1 • 10-14 


1.1 • 10-13 


1.7T 



Table 4: Comparison (absolute error, error estimation and computation time) of 
the algorithm of 3 (CC) and the present method (A) in the case of the integral 
/^^ log^(1.0001 - x){x - 0.99)-Ma; (IT = 0.0018s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


3.5 • 10-13 


2.9 • 10-13 


LOT 


A 


4.7-10-13 


9.0 • 10-12 


1.6T 



Table 5: Comparison (absolute error, error estimation and computation time) of 
the algorithm of 3 (CC) and the present method (A) in the case of the integral 
j\ VI cos(44x)|3(a; + 0.6)-^dx (IT = 2.3s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


7.3 • 10-12 


5.0 • 10-12 


LOT 


A 


4.9 • 10-15 


3.0 • 10-13 


0.007T 



Table 6: Comparison (absolute error, error estimation and computation time) of 
the algorithm of [3] (CC) and the present method (A) in the case of the integral 
J^^ VI -x^ cos(100a;)(a: - 0.9)-^dx (IT = 3.7s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


6.7-10-13 


3.7-10-12 


LOT 


A 


2.0 • 10-15 


7.1 • 10-14 


0.0009T 



for the widest possible class of functions. If the function / is not analytic 
(Tables [5] and [6]) , then the presented algorithm is much more efficient than 
the one of [3] . The results given in Tables [7] and [8] show that the problem 
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Table 7: Comparison (absolute error, error estimation and computation time) of 
the algorithm of [3] (CC) and the present method {A) in the case of the integral 
ill e'^ix - 0.9999999)-ida; (IT = 0.0003s). 



the algorithm 


absolute error 


error estimation 


computation time 


CC 


1.4 • 10-9 


1.4-10-15 


LOT 


A 


1.4 • 10-9 


3.0 • 10-9 


5.7T 



Table 8: Comparison (absolute error, error estimation and computation time) of 
the algorithm of j3] (CC) and the present method (A) in the case of the integral 
g-ioo(x+o.4)^ sin(e-i"^)(x + OAl)-^dx {IT = 0.003s). 



algorithm 


absolute error 


error estimation 


computation time 


CC 


1.3 • 10-13 


3.9 • 10-16 


LOT 


A 


1.2 • 10-13 


4.3 • 10-13 


0.3T 



discussed in Section [4. II is very important. When implementing the method 
of [3], we used only the error estimates presented there. As we can see, if the 
computations are performed in the fixed precision arithmetic, the additional 
error estimates similar to the ones given in (14. ID and (14. 2p should be a part 
of any algorithm for computing Cauchy principal value integrals. 

The experiments have shown that the very simple algorithm presented 
in this paper, based on the use of an adaptive quadrature, is very efficient 
and accurate. It can be applied to a very wide class of integrands and, as it 
was also proved theoretically, the influence of roundoff error on the accuracy 
of the result is very small. 
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